package revised_simplex;

import java.io.*;
import java.util.ArrayList;

import Jama.Matrix;

public class RevisedSimplex {
    public static void main(String[] args) {
        // 定义矩阵和向量
        ArrayList<ArrayList<Double>> A = new ArrayList<>();
        ArrayList<Double> b = new ArrayList<Double>();
        ArrayList<Double> c = new ArrayList<Double>();

        // 读入文件
        String filePath = "src/revised_simplex/sample.txt";
        try {
            FileReader fr = new FileReader(filePath);
            BufferedReader br = new BufferedReader(fr);
            String line;
            String[] temp;
            int i = 0;
            while ((line = br.readLine()) != null) {

                temp = line.split(" ");
                if (i == 0) {
                    // 填充c
                    for (int j = 0; j < temp.length - 1; j++) {
                        c.add(Double.valueOf(temp[j]));
                    }
                } else {
                    // 填充A和b
                    ArrayList<Double> tempList = new ArrayList<Double>();
                    for (int j = 0; j < temp.length; j++) {
                        tempList.add(Double.valueOf(temp[j]));
                        if (j == temp.length - 1) {
                            b.add(Double.valueOf(temp[j]));
                        }
                    }
                    A.add(tempList);
                }
                i++;

            }
            br.close();
        } catch (IOException e) {
        }
        // 实例化求解对象
        Solver solver = new Solver(A, b, c);
        solver.showModel();
        solver.solve();
    }
}

class Solver {
    ArrayList<ArrayList<Double>> A = new ArrayList<>();
    ArrayList<Double> b = new ArrayList<Double>();
    ArrayList<Double> c = new ArrayList<Double>();

    private ArrayList<Integer> x_b = new ArrayList<Integer>(); // 基变量
    private Matrix B_inverse; // 存储B逆
    private Matrix E; // 存储E，E左乘B逆得到新的B逆
    private ArrayList<Matrix> P = new ArrayList<Matrix>(); // 存储每次迭代的每个变量的技术参数向量
    private Matrix bb; // 和b一样，由于读文件的代码复制之前的，所以冗余了，不想改了
    private Matrix b_new; // 存储每次迭代的b'=B逆b
    private Matrix c_bv; // 基变量价值系数向量
    private boolean isOptimal = false; // 标识符，是否达到最优解

    private final int c_length;
    private final int b_length;

    Solver(ArrayList<ArrayList<Double>> A, ArrayList<Double> b, ArrayList<Double> c) {
        this.A = A;
        this.b = b;
        this.c = c;
        c_length = c.size();
        b_length = b.size();
    }

    void showModel() {
        // 输出目标函数
        String obj = "Max z = ";
        for (int i = 0; i < c_length; i++) {
            obj += c.get(i).toString() + "x" + Integer.toString(i);
            if (i != c_length - 1) {
                obj += " + ";
            }
        }
        System.out.println(obj);
        System.out.println("Subject to:");

        // 输出约束
        for (int i = 0; i < b_length; i++) {
            String con = "";
            for (int j = 0; j < c_length; j++) {
                con += Double.toString(A.get(i).get(j)) + "x" + Integer.toString(j);
                if (j != c_length - 1) {
                    con += " + ";
                }
            }
            con += " = " + Double.toString(b.get(i));
            System.out.println(con);
        }
    }

    private int getXin() {
        int xin = 0;

        // 计算c_bv
        for (int i = 0; i < b_length; i++) {
            c_bv.set(0, i, c.get(x_b.get(i)));
        }

        // 计算检验数
        double cmax = -Double.MAX_VALUE;
        for (int i = 0; i < c_length; i++) {
            // 更新变量技术参数向量
            Matrix p = new Matrix(b_length, 1);
            for (int j = 0; j < b_length; j++) {
                p.set(j, 0, A.get(j).get(i));
            }
            P.set(i, B_inverse.times(p));
            // 计算检验数
            double temp = c.get(i) - c_bv.times(P.get(i)).det();
            if (cmax < temp) {
                cmax = temp;
                xin = i;
            }
        }
        if (cmax <= 0) {
            isOptimal = true;
        }
        return xin;
    }

    private int getXout(int xin) {
        int xout = 0;
        double theta = Double.MAX_VALUE;
        // 更新b_new
        b_new = B_inverse.times(bb);
        // 计算theta
        for (int i = 0; i < b_length; i++) {
            double temp = b_new.get(i, 0) / P.get(xin).get(i, 0);
            if (temp < 0) {
                continue;
            }
            if (temp < theta) {
                theta = temp;
                xout = x_b.get(i);
            }
        }
        return xout;
    }

    private void updateB_inverse(int xin, int xout) {
        Matrix p = P.get(xin);
        // 构造E
        int index = x_b.indexOf(xout);
        for (int i = 0; i < b_length; i++) {
            Matrix temp = new Matrix(b_length, 1);
            if (i == index) {
                for (int j = 0; j < b_length; j++) {
                    if (j == index) {
                        temp.set(j, 0, 1 / p.get(j, 0));
                    } else {
                        temp.set(j, 0, -p.get(j, 0) / p.get(index, 0));
                    }
                }
            } else {
                for (int j = 0; j < b_length; j++) {
                    if (j == i) {
                        temp.set(j, 0, 1);
                    } else {
                        temp.set(j, 0, 0);
                    }
                }
            }
            E.setMatrix(0, b_length - 1, i, i, temp);
        }
        // 更新B_inverse
        B_inverse = E.times(B_inverse);
    }

    void solve() {
        // 确定初始基变量
        for (int i = 0; i < b_length; i++) {
            x_b.add(c_length + i - b_length);
        }

        // （假装）计算初始B逆
        double[][] B = new double[b_length][b_length];
        for (int i = 0; i < b_length; i++) {
            double[] temp = new double[b_length];
            for (int j = 0; j < b_length; j++) {
                temp[j] = A.get(i).get(x_b.get(j));
            }
            B[i] = temp;
        }
        B_inverse = new Matrix(B).inverse();

        // 生成初始E（单位矩阵）
        E = new Matrix(b_length, b_length);
        for (int i = 0; i < b_length; i++) {
            for (int j = 0; j < b_length; j++) {
                if (i == j) {
                    E.set(i, j, 1);
                } else {
                    E.set(i, j, 0);
                }
            }
        }

        // b_new初始等于b
        b_new = new Matrix(b_length, 1);
        for (int i = 0; i < b_length; i++) {
            b_new.set(i, 0, b.get(i));
        }
        bb = b_new.copy();

        // 填充c_bv
        c_bv = new Matrix(1, b_length);
        for (int i = 0; i < b_length; i++) {
            c_bv.set(0, i, c.get(x_b.get(i)));
        }

        // 填充P
        for (int i = 0; i < c_length; i++) {
            Matrix p = new Matrix(b_length, 1);
            for (int j = 0; j < b_length; j++) {
                p.set(j, 0, A.get(j).get(i));
            }
            P.add(B_inverse.times(p));
        }

        // 开始计算
        while (true) {
            B_inverse.print(3, 1);
            System.out.println("基变量为：" + x_b.toString());
            // 开始计算
            int xin = getXin();
            // 输出当前解
            System.out.println("Z = " + Double.toString(c_bv.times(b_new).det()));
            // 判断是否得到最优解
            if (isOptimal) {
                System.out.println("最优解为：");
                for (int i = 0; i < c_length; i++) {
                    System.out.println("x" + Integer.toString(i) + " = "
                            + (x_b.indexOf(i) != -1 ? Double.toString(b_new.get(x_b.indexOf(i), 0)) : "0"));
                }
                break;
            }
            int xout = getXout(xin);
            updateB_inverse(xin, xout);
            x_b.set(x_b.indexOf(xout), xin);
            System.out.println("==========================");
        }
    }
}